# level 1
function n_shapes0(ξ::Float64)
    n1 = ξ*(ξ - 1.0)/2.0; n2 = 1.0 - ξ*ξ; n3 = ξ*(ξ + 1.0)/2.0
    return [n1, n2, n3]
end
function n_shapes1(ξ::Float64)
    n1 = (2.0*ξ - 1.0)/2.0; n2 = -2.0*ξ; n3 = (2.0*ξ + 1.0)/2.0
    return [n1, n2, n3]
end

# level 2
function n_shapes(ξ::Float64, order::Int64)
    if order == 0
        return n_shapes0(ξ)
    elseif order == 1
        return n_shapes1(ξ)
    end
end

# level 3
function N_uc0(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[1,11,21]] = n_shapes(ξ,order); return f*r
end

function N_θc(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[2,12,22]] = n_shapes(ξ,order); return f*r
end

function N_wc(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[3,13,23]] = n_shapes(ξ,order); return f*r
end

function N_us0(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[4,14,24]] = n_shapes(ξ,order); return f*r
end

function N_us1(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[5,15,25]] = n_shapes(ξ,order); return f*r
end

function N_us2(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[6,16,26]] = n_shapes(ξ,order); return f*r
end

function N_us3(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[7,17,27]] = n_shapes(ξ,order); return f*r
end

function N_ws0(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[8,18,28]] = n_shapes(ξ,order); return f*r
end

function N_ws1(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[9,19,29]] = n_shapes(ξ,order); return f*r
end

function N_ws2(ξ::Float64, order::Int64, xi::Float64, xj::Float64)
   invj = 2.0/(xj - xi); f = invj^order
   r = zeros(1,30); r[[10,20,30]] = n_shapes(ξ,order); return f*r
end

# level 4
function N_Uc(ξ::Float64, order::Int64, z::Float64, xi::Float64, xj::Float64)
   r = N_uc0(ξ,order,xi,xj) - z*N_θc(ξ,order,xi,xj)
   return r
end
function N_Wc(ξ::Float64, order::Int64, z::Float64, xi::Float64, xj::Float64)
   r = N_wc(ξ,order,xi,xj)
   return r
end
function N_Us(ξ::Float64, order::Int64, z::Float64, xi::Float64, xj::Float64)
   r = N_us0(ξ,order,xi,xj) + z*N_us1(ξ,order,xi,xj) + z*z*N_us2(ξ,order,xi,xj) + z*z*z*N_us3(ξ,order,xi,xj)
   return r
end
function N_Ws(ξ::Float64, order::Int64, z::Float64, xi::Float64, xj::Float64)
   r = N_ws0(ξ,order,xi,xj) + z*N_ws1(ξ,order,xi,xj) + z*z*N_ws2(ξ,order,xi,xj)
   return r
end

# level 5, defining shape functions for strains
function N_ϵcx(ξ::Float64, z::Float64, xi::Float64, xj::Float64)
   return N_Uc(ξ, 1, z, xi, xj)
end
function N_γc(ξ::Float64, xi::Float64, xj::Float64)
   r = N_wc(ξ, 1, xi, xj) - N_θc(ξ,0, xi, xj)
   return r
end
function N_ϵsx(ξ::Float64, z::Float64, xi::Float64, xj::Float64)
   return N_Us(ξ, 1, z, xi, xj)
end
function N_ϵsz(ξ::Float64, z::Float64, xi::Float64, xj::Float64)
   r = N_ws1(ξ,0,xi,xj) + 2.0*z*N_ws2(ξ,0,xi,xj); return r
end
function N_γs(ξ::Float64, z::Float64, xi::Float64, xj::Float64)
   p1 = N_us1(ξ,0,xi,xj) + 2.0*z*N_us2(ξ,0,xi,xj) + 3.0*z*z*N_us3(ξ,0,xi,xj)
   p2 = N_ws0(ξ,1,xi,xj) + z*N_ws1(ξ,1,xi,xj) + z*z*N_ws2(ξ,1,xi,xj)
   r = p1 + p2
   return r
end
function N_ucs(ξ::Float64, xi::Float64, xj::Float64,h1::Float64, h2::Float64)
   r = N_Uc(ξ,0,-h1,xi,xj) - N_Us(ξ,0,h2,xi,xj); return r
end
function N_wcs(ξ::Float64, xi::Float64, xj::Float64,h1::Float64, h2::Float64)
   r = N_Wc(ξ,0,-h1,xi,xj) - N_Ws(ξ,0,h2,xi,xj); return r
end

# level 6, forming strain matrices B
function B_matrix_c(ξ::Float64, z::Float64, xi::Float64, xj::Float64)
   return [N_ϵcx(ξ, z, xi, xj); N_γc(ξ, xi, xj)]
end

function B_matrix_s(ξ::Float64, z::Float64, xi::Float64, xj::Float64)
   return [N_ϵsx(ξ, z, xi, xj); N_ϵsz(ξ, z, xi, xj); N_γs(ξ, z, xi, xj)]
end
